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Abstract. In this study, we propose an alternative way to simulate particle 
suspensions using the lattice Boltzmann method. The main idea is to impose the 
non-slip boundary condition in the lattice sites located on the particle boundaries. 
The focus on the lattice sites, instead of the links between them, as done in the 
more used methods, represents a great simplification in the algorithm. A fully 
description of the method will be presented, in addition to simulations comparing 
the proposed method with other methods and, also, with experimental results. 
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1. Introduction 

Particulate flows are found in many industrial processes [H [2 |3l [4] , and have been 
subjected to considerable scientific investigation. Recently, computer simulations have 
become an effective tool in these studies and several methods have been applied, 
such as finite element method [5] , Lagrange- multipliers [H |3] , direct forcing method 
[6], lattice Boltzmann methods [Il[l[7l[a[i[10l[IIl[Il|ia[ll[l3[16], solving the 
Stokes equation near the particles (Physalis) [17] and combining two or more methods 

In this study, we propose an alternative way to simulate suspensions using the lattice 
Boltzmann method. The main idea is to impose the non-slip boundary condition in 
the lattice sites representing the particle boundaries. The focus on lattice sites, instead 
of the link between them, as done in the more used methods [11 , represents a great 
simplification in the algorithm. Similar approaches, focusing on the lattice sites were 
already tried p. , although without popularity due, probably, to difficulties several in 
the implementation. 

Its important to mention that the simplifications we propose can reduce the accuracy 
in describing the details of the flow near the particles, and these details can or cannot 
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Figure 1. The D2Q9 lattice, used in two-dimensional simulations. The arrows 
represent the distributions fi(x,t), and circles represent the lattice sites. 

be important, depending on the problem one wants to simulate. In the simulations we 
have done until now, in despite of the simplifications, the results we obtain are similar 
to the results obtained by other methods (see section |3|. 

2. The model 

In this section we introduce the lattice Boltzmann method and present the model 
describing the particle-fiuid interaction. Interactions among particles and between a 
particle and a solid surface will also be discussed in this section. 

2.1. The lattice Boltzmann model 

The lattice-Boltzmann method (LBM) is based on the discretization of the 
Boltzmann's mesoscopic equation[20l IHJ [22], usually with the BGK approach for 
the collision operator (for a comprehensive review see [23 ). In the LBM scale, the 
system is described using a single particle distribution function, /^(ic,^), representing 
the number of particles with velocity q at the site x and time where z = 0, 6. The 
particles are restricted to a discrete lattice, in a manner that each group of particles can 
move only in a finite number b of directions and with a limited number of velocities (see 
Fig. [T]). Therefore, physical and velocity space are discretized. The local macroscopic 
properties such as total mass (the particle mass, m, is assumed unitary), p(ic,t), and 
total momentum, p{x^t)u{x,t)^ can be obtained from the distribution function in the 
following way: 

p{x,t) = j2fiix,t), (1) 

i 

p{x, t)u{x, t) = Y^ fi{x, t)ci. (2) 

i 

The Lattice Boltzmann equation, that is, the discrete version of the Boltzmann 
equation with the BGK collision, operator is written 

Mx + Stcut + 5t) -Mx,t) = -^-^ [Mx,t) - f:\p,u)] (3) 
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Figure 2. Boundary sites characterizing a particle. The boundary sites (BS) are 
depicted in dark gray and the internal sites (IS), in dark. The particle's boundary 
are regard to be between the IS and BS sites. 

where, 6t is the time step and fi^{p, u) is a polynomial approximation of the Maxwell- 
Boltzmann equilibrium distribution [24l |20l [2TJ |22] , a function of the local variables 
p{x^t) and u{x^t). It can be shown, through a Chapman-Enskog analysis, that this 
system macroscopically will evolve according to the Navier- Stokes equations with a 
kinematic viscosity given by 

p = cUt- 1/2) , (4) 

where Cs is the sound velocity, a constant depending on the set of velocities 

2.2. Particle- fluid interaction 

The basic idea of the method is that the fluid in contact with a solid surface must 
acquire the velocity of this surface, considering the non-slip condition. Keeping this 
in mind, a set of boundary sites can be used to describe the particle. This approach 
is similar to the one presented in and [11], although we focus in the lattice sites 
and not in the links between them. We denote boundary sites (BS) the particle sites 
in contact with fluid, and internal sites (IS) the particle sites in contact with the 
boundary sites (see Figj2]). The particle's boundary is regarded to be halfway between 
the IS and the BS sites. A particle will be represented by its central position Xp^ a 
radius r^, a mass nip and a moment of inertia Ip. Its velocity will be denoted i/p, and 
the angular velocity ujp. Only spherical particles will be considered. 

The presence of a particle will be represented only by its effects on the fluid in 
contact with the particle, altering the dynamics of the sites BS and IS. The position of 
these sites will be denoted hy Xb and cc/, respectively, and its integer counterpart (or 
the nearest integer of each component of the vectors) by [xb] and [x/], respectively. 
Velocities and densities in the sites [xb] and [xj] will be denoted ub^ uj and p^, 
Pi. We emphasize that all sites will be updated by the usual collision/propagation 
steps, but the IS and BS sites, in addition to these steps, will be submitted to a set 
of substeps included between the collision and the propagation steps. These substeps 
are presented in the sequence. 
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• Particle-Buid momentum transfer 

This substep is executed only in the BS sites. As was previously mentioned, the 
velocity in the BS sites must acquire the particle's velocity. This velocity, considering 
the angular velocity, is 

u^^ = Up^ {xb - Xp) X ujp. (5) 

To impose this velocity in the fluid, that is in lattice sites, the equilibrium 
distribution is employed, setting fi {[xb] ,t) = f^^{p^u'^). 

Naming ps = PbUb the linear momentum of a BS site, we compute 

Aps = Pb{ub -ub), (6) 

and 

AZb = Pb{xb - Xp) X {u'b - ub), (7) 

where /SI b represents the change in the angular momentum caused by the change in 
the linear momentum Ap^. 

Finishing this step we have the total momentum exchanged, APg and ALb- 



APb = J2^Pb (8) 

BS 

ALb = J2^^b (9) 

BS 

• Particle^ s acceleration 

This substep accounts for the change in particles velocity caused by fluid. 
According to the Newton's third law of motion, we simply compute 

u'^ = Up- APb/ttIp, (10) 

uj'p=i^p- ALb/Ip, (11) 

u'j = Up^ {xi - Xp) X ujp. (12) 



Updating of particle^ s position 



Due to the spherical symmetry it is not necessary to take into account the rotation 
of the particles. The position is updated by doing: 

Xp = Xp-\- StUp. (13) 

The positions of boundary and internal sites, also, must be updated: 

Xb = xb + ^tu'p] (14) 

x'j = xi -\- 6tUp. (15) 

Clearly, this procedure includes errors of order (^|), more precise procedures could 
be employed. Nevertheless, there is a lot of imprecision in describing the shape of the 
particle, therefore it is not necessary to be so precise in updating particle's position. 

• Velocity of the internal sites 
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To impose this velocity the equihbrium distribution is employed again, setting 
fi{[x'A,t) = f:\p,u',). 

Some comments are necessary to clarity the physics behind these sub-steps. Con- 
sider a particule that is found with a velocity Up different of the fluid velocity around 
it, in the beginning of the particle-fluid interaction step. This situation never occurs 
in the continuum but this is possible in discrete models due to the discretization of 
time. Particle and fluid must then exchange linear and angular momentum until the 
fluid around the particle acquires the velocity of the particle surface. In this pro- 
cess the particle velocity also changes due to the fluid reaction on the particle and in 
accordance with Newton's third law of motion. Action and reaction happen simul- 
taneously in the physical process, but in the discrete case this occurs in steps 1 and 
2. In the first step the particle transfers linear and angular momentum to the fluid. 
In consequence the fluid in contact with the particle surface particle accelerates until 
have acquired the same particle velocity. As it was earlier mentioned two hypotheses 
are used for this step: non-slipping on the particle surface and the hypotheses that 
the fluid-particle equilibration takes a time interval smaller than the time step used 
in the simulation. This enables the use of an equilibrium distribution to impose the 
velocities. With the fluid velocity altered, we then calculate the change in the linear 
and angular momentum of the fluid due to this acceleration. These variations must 
be the same as the corresponding changes for the particle, in accordance with the 
Newton's third law. We then proceed to step 2 and recalculate the particle velocity. 
The velocity to be imposed in sites IS is also calculated in accordance with the new 
particle velocity. 

In a third step, the position of the particle and sites IS and BS are changed. The 
spherical symmetry of the particles eases this step since only the translation velocities 
are taken into account in the calculation of the new positions. It must be stressed 
that although sites BS belong to the fluid phase, they have their position changed in 
accordance with the particle velocity. Indeed, these sites are in the particle boundary 
and must follow its displacement. In other words the displacement of the boundary 
sites is independent of the fluid particles that are occupying these sites, in a given 
instant. 

Finally, in a fourth step, we impose the velocity Uj to the sites IS. This step 
is important because during the propagation step the information in these sites 
is transferred to the adjacent sites BS. Therefore, the sites BS will have their 
velocity composed by the particle and fluid velocities and will define a new change of 
momentum in the following time step. 

2. 3. Interaction between particles 

Especially in simulations involving a great number of particles, as the simulation 
presented in subsection |3.4[ the interaction between particles must carefully be ad- 
dressed. There are many techniques to treat these interactions, possibly the most 
popular approach is to introduce a repulsive force between particles when the gap be- 
tween becomes smaller than a given threshold [14 . We chose another approach in this 
work, treating the collisions as occurring between completely rigid particles (hard-core 
collisions). In this way it is not necessary to set the parameters of repulsive forces. 
We simply impose (see Figj3| 
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Figure 3. Collision between two particles, A and B. The arrows indicate the 
normal and tangential velocities before and after the collision. The velocities are 
represented in the center of mass reference frame. 



'^'au = -'^An^ (16) 
' = -VBn. (17) 

VAU (18) 
'^Bt = '^Bt^ (19) 

where all velocities are represented in the center of mass reference frame. This 
imposition must be done before updating particles' positions, that is, it can be 
considered as part of the particle's acceleration substep. 

It's necessary to underline that the way we chose to introduce the particle-particle 
interaction is not new, nor regarded as part of the method. Moreover, this choice 
was done considering solely the easeness in the implementation. In accordance with 
the Reynolds number, volume fraction of the particles and the involved force, more 
elaborated techniques, employing lubrication [8 or spring forces [25 may reveal to be 
necessary. 



3. Validation 



Four cases are presented in order to validate the model. The first case simulated 
was the flow around a massive two-dimensional particle. The particle's mass was 
chosen in such way that it doesn't move, therefore it's possible to compare the flow 
around the two dimensional particle and the flow around a solid cylinder, simulated 
using bounce-back. The drag coefficient for a massive particle was also computed 
considering several Reynolds number and compared with drag coefficients obtained 
by other methods. The second case simulated was the flow of a neutrally buoyant 
two-dimensional particle in a shear flow, the results are compared with the results 
presented by Feng and Michaelides[14]. A simulation of a sphere settling in a closed 
box was done in order to validate the model in a three-dimensional simulation. The 
results obtained was compared with the ones presented by ten Gate et al.[26^. Finally, 
it was simulated the sedimentation of 504 two-dimensional particles in an enclosure. 
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Figure 4. Geometry used to simulate the flow in a channel with an obstacle. 



3.1. Flow around a massive two dimensional particle 

The first geometry used to simulate tiie flow around a massive particle is depicted 
in Fig. |4] The simulations of a flow around a solid cylinder was carried out using 
bounce-back boundary conditions and the simulation of flow around a massive particle 
was done using the method proposed in this work, that is, using the equilibrium 
distribution to impose the velocities in the boundaries of the particle. The results of 
both simulations are presented in Fig. [5) To emphasize the deviations we plot in Fig. 
[6] the magnitude of the difference between the velocities. In the enclosed regions of 
Fig. [6] the deviations are of order 10 ^ (the velocities in the simulations varying from 
zero to 0.075). Is important to notice that, in despite of deviations that can appear 
near the solid surface as result of applying different boundary conditions, the overall 
flow behavior is the same. Off course, depending on the applications ones intends to 
focus, these differences must be taken into account. 

Simulations to obtain drag coefficients were also carried out in the geometry 
presented in Fig. [7j Velocities varying from 0.00417 to 0.1667 in the x direction were 
imposed at the boundaries (left, right, top and bottom), resulting in the Reynolds 
number varying from Re = 1 to Re = 40. The results obtained are shown in Fig. 
[sj beside the results published by Rajani et al.[27 and Silva et al. [28 . Although 
the drag coefficients computed in this work were systematically lower than the results 
obtained by other methods, the importance of these errors is dependent on what we 
want to describe in a given problem, as it can be seen in the simulations presented in 
the sequence. 



3.2. Neutrally buoyant two dimensional particle in a shear flow 

The motion of a neutrally buoyant two-dimensional particle moving in viscous fluid 
was already simulated using LBM [29 ,[30 ,[14 as well as using finite element method 
[5] , and was chosen as one of the validations of the present model. The geometry 
of the problem is described in Fig. [oj where /7^(;/2 and —Uw/2 are the velocities 
imposed. Periodic boundary conditions are imposed in the left and right boundaries. 
The relaxation time was set r = 0.6, which implies a kinematic viscosity v = 1/30, 
in lattice units. The parameters chosen are the same used in Feng & Michaelides 
paper [14 , in order to ease the comparison. The velocity equals Uyj/2 = 1/120, 
therefore, the shear rate for the flow is 7 = U^/H = 1/120, and a dimensionless 



time, t*, is deflned = tj'^rs/i'^ where is the particle radius. Fig. 10, shows the 



migration of the particle, initially at the position y = i^/4, toward the center. The 
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Figure 6. The difference in the velocities considering the simulations using a 
solid cylinder and a massive particle. In the enclosed region the errors are of 
order 0(10""^), outside the errors are smaller than this. 
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Figure 7. The geometry used to compute the drag coefficient of the 2D particle. 
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Figure 8. Drag coefficients for Reynolds number varying from one to forty. The 
results computed using the proposed model are compared with the work by Rajani 
et al. (27 (Reynolds from one to five) and Silva et al. 28 (Reynolds from ten to 
forty). 
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Figure 9. Schematic depiction of the problem simulated: a neutrally buoyant 
two dimensional particle in a shear flow 
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Figure 10. A neutrally buoyant particle in a shear flow - comparison of the 
results obtained using the proposed model for fluid-solid interaction and previous 
models. The y-axis shows the particle's position with respect to the lateral of 
channel divided it's width H. 



agreement between the results using the proposed model and the previous models are 
quite good. 



3.3. A sphere settling in a closed box 

In this subsection the trajectory and velocity of a sphere settling in a closed box 
is simulated using the model here prosed and compared with experimental results 
obtained by Gate et al. The settling sphere has a diameter oi D = 15mm 

and density p = 1120kg / m^ . The container dimensions are depth x width x height = 
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Figure 11. Schematic depiction of the geometry used in the experiments (left) 
and in the simulations (right). 



100 X 100 X 160mm (see Fig. 11). Four cases were simulated considering the different 



densities and viscosities of the fluid in which the sphere will settle, corresponding to a 
Reynolds number varying from Re = 1.5 to Re = 32.2. The fluid characteristics and 
the parameters used in the simulations are shown in Table [T] 



Table 1. Fluid properties in the experiment and parameters used in simulations. 





Pf[kg/m^] 


/j.fllO-'-^Ns/m^] 


r 


St[lO--'s] 


Case El 


970 


373 


0.9 


0.347 


Case E2 


965 


212 


0.9 


0.607 


Case E3 


962 


113 


0.8 


0.851 


Case E4 


960 


58 


0.9 


2.207 



Figure [12] and [13] present comparisons between simulation and experiment for the 
trajectories and velocities of the settling particle, respectively. 

3-4- 504 particles settling in a closed box 

The problem of a large number of particles settling in a closed 2D box was already 
simulated by other methods 0, [TJ- All the parameters were chosen in order to 
compare with the previous works. That is, 504 circular particles with diameter 
D = 0.625cm settling in box having 10~^m width and 10~^m height. The fluid 
and particle densities are Pf = lOOOkg/m^ and pp = 1010/c^/m^, respectively, and 
the fluid kinematic viscosity is Uf = 10~^m^/5. Representing the box by 512 x 512 
sites and using a relaxation time r = 0.9915 we will have a time step of 0.00025s. 
The process of sedimentation simulated from the initial state to 245 is presented in 
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Figure 12. The position of the particle setthng in a closed box. The simulated 
results are compared with the experimental results from Gate el al. (2002). The 
y-axis presents particle's position H (see Fig. divided by it's diameter D. The 
X-axis shows time in seconds. 
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Figure 13. Comparison between simulation and the experimental results from 
Gate el al. (2002). The y-axis shows the velocity of the settling particle in meters 
per second and the x-axis indicates time in seconds. 
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Figure 14. Configurations of the settling particles. 



t=2s 




t=3s t=4s 



Figure 15. Configurations of the settling particles. 
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Figure 16. Configurations of the settling particles. 



t=24s 



Fig. [T4j Fig. [15] and Fig. [T6j The figures show the development of a Rayleigh- 
Taylor instabihty and are, quahtatively, in agreement with the previous works. The 
differences between the three simulations (the one presented here and the simulations 
of refs. ^ and [14 ) are, possibly, a result of the differences in treating the collisions 
between particles and differences arising from compressibility effects that are present 
in the lattice Boltzmann methods. 
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4. Conclusion 

In this study an alternative and simpler way to simulate particle-fluid interactions is 
proposed. The lattice Boltzmann method is employed to simulate the fluid flow and 
the particles are simulated using the Newton's law. The coupling is made applying 
the equilibrium distribution function to assure the non-slip condition near the solid 
surfaces. Several simulations are presented showing that the method can simulate 
particle-fluid interactions with a precision comparable with other methods. 
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